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We introduce novel algorithms for the quantum simulation of molecular systems which are asymp¬ 
totically more efficient than those based on the Trotter-Suzuki decomposition. We present the first 
application of a recently developed technique for simulating Hamiltonian evolution using a truncated 
Taylor series to obtain logarithmic scaling with the inverse of the desired precision, an exponential 
improvement over all prior methods. The two algorithms developed in this work rely on a second 
quantized encoding of the wavefunction in which the state of an N spin-orbital system is encoded in 
0{N) qubits. Our first algorithm requires at most 0{N^t) gates. Our second algorithm involves on- 
the-fly computation of molecular integrals, in a way that is exponentially more precise than classical 
sampling methods, by using the truncated Taylor series simulation technique. Our second algorithm 
has the lowest gate count of any approach to second quantized quantum chemistry simulation in the 
literature, scaling as 0{N^t). The approaches presented here are readily applicable to a wide class 
of fermionic models, many of which are defined by simplified versions of the chemistry Hamiltonian. 


I. INTRODUCTION 

As small, fault-tolerant quantum computers come in¬ 
creasingly close to viability [1-4] there has been substan¬ 
tial renewed interest in quantum simulating chemistry 
due to the low qubit requirements and industrial impor¬ 
tance of the electronic structure problem. A recent se¬ 
ries of papers tried to estimate the resources required 
to quantum simulate a small but classically intractable 
molecule [5-9]. Although qubit requirements seem mod¬ 
est, initial predictions of the time required were daunting. 
Using arbitrarily high-order Trotter formulas, the tight¬ 
est known bound on the gate count of the second quan¬ 
tized, Trotter-based quantum simulation of chemistry is 
0{N^t/e°^^'^) [10, 11]^, where e is the precision required 
and N is the number of spin-orbitals. However, using 
significantly more practical Trotter decompositions, the 
best known gate complexity for this quantum algorithm 
isOiN^y/t^e) [ 6 ]. 

Fortunately, numerics indicated that the average 
circuit depth for real molecules may be closer to 
d{N^y/¥Je) [7], or s/WJl) when only trying 

to simulate ground states, where Z is the largest nuclear 
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charge for the molecule [9]. While this improved scal¬ 
ing restores hope that fault-tolerant devices will have an 
impact on some classically intractable chemistry prob¬ 
lems, the Trotter-based quantum simulation of large 
(e.g. N > 500) molecules still seems prohibitively costly 
[9, 12, 13]. This limitation would preclude simulations 
of many important molecular systems, such as those in¬ 
volved in biological nitrogen fixation and high-T), super¬ 
conductivity [12, 13]. 

The canonical quantum algorithm for quantum chem¬ 
istry, based on the Trotter-Suzuki decomposition which 
was first applied for universal quantum simulation in 
[14, 15], was introduced nearly one decade ago [16]. This 
approach was later refined for implementation with a set 
of universal quantum gates in [17]. With the exception 
of the adiabatic algorithm described in [18] and a clas¬ 
sical variational optimization strategy making use of a 
quantum wavefunction ansatz described in [19], all prior 
quantum algorithms for chemistry have been based on 
Trotterization [20-27] . 

Trotter-Suzuki approaches were also applied to simu¬ 
lation of evolution under sparse Hamiltonians with the 
entries given by an oracle [28, 29]. A related problem is 
the simulation of continuous query algorithms; in 2009, 
Cleve et al. showed how to achieve such simulation with 
exponentially fewer discrete queries than Trotterization 
in terms of 1/e [30] . The algorithm of [30] still required a 
number of ancilla qubits that scaled polynomially in 1/e, 
but this limitation was overcome in [31] which demon¬ 
strated that the ancilla register in [30] could be com¬ 
pressed into exponentially fewer qubits. In [32, 33], Berry 
et al. combined the results of [28-31] to show exponen¬ 
tially more precise sparse Hamiltonian simulation tech¬ 
niques. A major contribution of [32] was to use oblivi- 
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ous amplitude amplification to make the algorithm from 
[30, 31] deterministic, whereas prior versions had relied 
on probabilistic measurement of ancilla qubits. An im¬ 
provement introduced in [33] was to show how to sim¬ 
ulate arbitrary Hamiltonians using queries that are not 
self-inverse (a requirement of the procedure in [32]). We 
focus on the methodology of [33] which is relatively self- 
contained. 

The algorithm of [33] approximates the propagator us¬ 
ing a Taylor series expansion rather than the Trotter- 
Suzuki decomposition. By dividing the desired evolution 
into a number of simulation segments proportional to the 
Hamiltonian norm, one can truncate the Taylor series at 
an order which scales logarithmically in the inverse of the 
desired precision [33]. The truncated Taylor series must 
be expressed as a weighted sum of unitary operators. To 
simulate the action of this operator, one first initializes 
the system along with an ancilla register that indexes 
all terms in the Taylor series sum. The ancilla register is 
then put in a superposition state with amplitudes propor¬ 
tional to the coefficients of terms in the Taylor series sum. 
Next, an operator is applied to the system which coher¬ 
ently executes a single term in the Taylor series sum that 
is selected according to the ancilla register in superposi¬ 
tion. Finally, by applying the transpose of the procedure 
which prepares the ancilla register, one probabilistically 
simulates evolution under the propagator. The algorithm 
is made deterministic using an oblivious amplitude am¬ 
plification procedure from [32]. 

This is the first paper of a two-paper series which ap¬ 
plies the algorithm of [33] to quantum chemistry simu¬ 
lation. The algorithms discussed in this paper employ 
a second quantized encoding of the Hamiltonian, where 
we dynamically perform the Jordan-Wigner transforma¬ 
tion [34, 35] on the quantum computer. In the second 
paper of this series, we use a compressed, first quantized 
encoding of the wavefunction which requires a number 
of qubits that scales almost linearly with the number of 
electrons [36]. 

In the present paper we develop two new algorithms 
for the application of the Hamiltonians terms, which we 
refer to as the “database” algorithm and the “on-the- 
fly” algorithm. In the database algorithm, the ancilla 
register’s superposition state is prepared with amplitudes 
from a precomputed classical database. In the on-the-fly 
algorithm, those amplitudes are computed and prepared 
on-the-fly, in a way that is exponentially more precise 
than classically possible. 


II. OVERVIEW OF RESULTS 

The simulation procedure described in [33] assumes the 
ability to represent the Hamiltonian as a weighted sum 
of unitaries which can be individually applied to a quan¬ 
tum state. Specifically, we must be able to express the 


simulation Hamiltonian as 

r 

7=1 

where the W-y are complex-valued scalars^, the are 
unitary operators and a mechanism is available for selec¬ 
tively applying the H^. Using the Jordan-Wigner trans¬ 
formation [34, 35] or the Bravyi-Kitaev transformation 
[37—39], the second quantized molecular Hamiltonian can 
be mapped to a sum of T € 0{N^) local Hamiltonians. 
Since these local Hamiltonians are each a tensor product 
of Pauli operators multiplied by some coefficient, they 
automatically satisfy the form of Eq. (1). 

We will need a circuit referred to in [33] as SELECT(i7) 
which is queried within the algorithm such that 

SELECT (H) |7) IV') = l7) IV’) ■ (2) 

One could construct SELECT(i7) by storing all the Pauli 
strings in a database. However, accessing this data would 
have time complexity of at least H(r). Instead, we com¬ 
pute and apply the Pauli strings using 0{N) gates (which 
can be parallelized to 0{1) circuit depth) by dynami¬ 
cally performing the Jordan-Wigner transformation on 
the quantum hardware. 

The algorithm of [33] also requires an operator that we 
refer to as prepare (W) which applies the mapping 

prepare(W)| 0)®‘°®^ = 7i^yPu;|7) (3) 

^ 7=1 

where 

r 

A = ^|VF^|, A&0{N^) (4) 

7=1 

is a normalization factor that will turn out to have signif¬ 
icant ramifications for the algorithm complexity. In the 
first of two algorithms discussed in this paper, we im¬ 
plement prepare(IU) using a database via a sequence 
of totally controlled rotations at cost 0(r). Because our 
first approach uses a database to store classically precom¬ 
puted values of Wj in order to implement prepare(JU), 
we refer to the first algorithm as the “database” algo¬ 
rithm. 

While we suggest a different strategy in Section HI, 
a database could also be used to construct SELECT(i7). 


^ The convention of [33] requires that the Wj are real, non-negative 
scalars. This treatment remains general as arbitrary phases can 
be factored into the Hj. However, we break with that convention 
and allow the W-y to take arbitrary complex values. This is done 
for pedagogical purposes: so that we may separately describe 
computation of the H^y and the W-y for the chemistry Hamilto¬ 
nian. Consequentially, our Eq. (41) differs from the analogous 
equation in [‘^3] by a complex conjugate operator. 
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That is, a controlled operation is performed which ap¬ 
plies Hi if 7 = 1 , followed by a controlled operation 
which performs H 2 if 7 = 2, and so forth. This would 
result in a slightly higher gate count than PREPARe(IT), 
because each of the T controlled operations must act on 
O(logiV) qubits even if the Bravyi-Kitaev transforma¬ 
tion is used. Nevertheless, this might represent a simpler 
solution than our construction of SELECt( 77) for early 
experimental implementations. 

Our second algorithm involves modifications to the al¬ 
gorithm of [33] which allows us to avoid some of this 
overhead. We exploit the fact that the chemistry Hamil¬ 
tonian is easy to express as a special case of Eq. (1) in 
which the coefficients are defined by integrals such as 

H4y = J w-y (z) dz. (5) 

Because our approach involves computing integrals on- 
the-fly, we refer to the second algorithm as the “on-the- 
fly” algorithm. We begin by numerically approximating 
the integrals as finite Riemann sums such as 

V . ^, 

^ (zp) (6) 

^ P=i 

where Zp is a point in the integration domain at grid 
point p. Equation ( 6 ) represents a discretization of the 
integral in Eq. (5) using p, grid points where the domain 
of the integral, denoted as Z, has been truncated to have 
total volume V. This truncation is possible because the 
functions w^{z) can be chosen to decay exponentially for 
the molecular systems usually studied in chemistry. Note 
that this might not be true for other systems, such as 
conducting metals. 

Our algorithm is effectively able to numerically com¬ 
pute this integral with complexity logarithmic in the 
number of grid points. It might be thought that this 
is impossible, because methods of evaluating numeri¬ 
cal integrals on quantum computers normally only give 
a square-root speedup over classical Monte-Carlo algo¬ 
rithms [40] . The difference here is that we do not output 
the value of the integral. The value of the integral is only 
used to control the weight of a term in the Hamiltonian 
under which the state evolves. 

We construct a circuit which computes the values of 
w-y{zp) for the quantum chemistry Hamiltonian with 
0{N) gates. We call this circuit SAMPLe(i(;) and define 
it by its action, 

SAMPLE {w) I 7 ) Ip) JO)®“ = I 7 ) Ip) 1 (fp)), (7) 

where w-y{zp) is the binary representation of w-y{zp) using 
logM qubits. 

By expanding the Wj in Eq. (1) in terms of the easily 
computed Wy{z} as in Eq. ( 6 ), we are able to compute 
analogous amplitudes to those in Eq. (3) in an efficient 
fashion. Thus, we no longer need the database that char¬ 
acterizes that algorithm. State preparation where the 


state coefficients can be computed on the quantum com¬ 
puter is more efficient than when they are stored on, and 
accessed from, a database [41]. The worst-case complex¬ 
ity is the square root of the dimension (here it would be 
O(vTp)), whereas the database state preparation has 
complexity linear in the dimension (which is 0 (r) for 
Hfy). Here this would not be an improvement, as we 
have increased the dimension in the discretization of the 
integral. 

However, the worst-case complexity is only if the am¬ 
plitudes can take arbitrary values (as this would enable a 
search algorithm, where the square root of the dimension 
is optimal [42]). If the amplitudes differ only by phases, 
then complexity of the state preparation is logarithmic 
in the dimension. We therefore decompose each w^{z) 
into a sum of terms which differ only by a sign. Then, 
although the dimension is increased, the complexity of 
the state preparation is reduced. The decomposition is 
of the form 


M 

Wy (Z) « C X! ^ 

m—1 

where 

Mg 0 ^maxlu;^(z)l/C^ . (9) 

In turn, we can express the Hamiltonian as a sum of 
unitaries weighted by identical amplitudes which differ 
only by an easily computed sign. 


r M II 

H = — W-f^m (Zp) Hj. 

^ 7 = 1 771=1 p—1 


( 10 ) 


As discussed above, the state preparation needed can 
be performed much more efficiently because the ampli¬ 
tudes are now identical up to a phase. By making a single 
query to SAMPLe(i/;) and then performing phase-kickback 
we can implement the operator prepare(i(;) whose ac¬ 
tion is 


PREPARE (u>) [0)®'°®^^^ = I 


where \i) = [ 7 ) \m) \p), L G Q{rMp) and 
A = L— G 0 ( TV rnax (z)| 




^,1 


( 11 ) 

( 12 ) 


is a normalization factor that will turn out to have signif¬ 
icant ramifications for the algorithm complexity. Later, 
we will show that A G 0{N^) and that prepare(w) can 
be implemented with 0{N) gate count, the cost of a sin¬ 
gle query to SAMPLE(r(;). 

The database algorithm performs evolution under H 
for time t by making 0{At) queries to both SELECT(i7) 
and prepare(IE). Because prepare(IE) requires 0(r) 
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= 0{N'^) gates, the overall gate count of this ap¬ 
proach scales as 0{N'^M). To avoid the overhead from 
prepare(IT), our on-the-fly algorithm exploits a mod- 
ihed version of the truncated Taylor series algorithm 
which allows for the same evolution by making 0{Xt) 
queries to SELECT(i7) and prepare(w). As PREPARE(?ii) 
requires 0{N) gates, the gate count for our on-the-fly al¬ 
gorithm scales as 0{NXt). 

The paper is outlined as follows. In Section III we in¬ 
troduce the second quantized encoding of the wavefunc- 
tion and construct select (i?). In Section IV we review 
the procedure in [33] to demonstrate our database algo¬ 
rithm which uses SELECT(iJ) and prepare(IT) to per¬ 
form a quantum simulation which is exponentially more 
precise than Trotterization. In Section V we show that 
one can modify the procedure in [33] to allow for es¬ 
sentially the same result while simultaneously comput¬ 
ing the integrals on-the-fly, and show how to implement 
prepare(ui) so as to compute the integrals on-the-fly. 
In Section VI we bound the errors on the integrals by 
analyzing the integrands. In Section VII we discuss ap¬ 
plications of these results and future research directions. 


III. THE HAMILTONIAN ORACLE 


The molecular electronic structure Hamiltonian de¬ 
scribes electrons interacting in a fixed nuclear potential. 
Using atomic units in which the electron mass, electron 
charge. Coulomb’s constant and h are unity we may write 
the electronic Hamiltonian as 




\R^ - Tj 


y 

It- — ■ 


(13) 


i,j>i 


where Ri are the nuclei coordinates, Zi are the nu¬ 
clear charges, and fi are the electron coordinates [43]. 
We represent the system in a basis of N single-particle 
spin-orbital functions usually obtained as the solution 
to a classical mean-held treatment such as Hartree-Fock 
[43]. Throughout this paper, (pi(rj) denotes the spin- 
orbital occupied by the electron which is parameter¬ 
ized in terms of spatial degrees of freedom fj. 

In second quantization, antisymmetry is enforced by 
the operators whereas in hrst quantization antisymmetry 
is explicitly in the wavefunction. The second quantized 
representation of Eq. (13) is 

H = 'y hija\aj + ^ X! hijkia\a]akai (14) 

ij ijk£ 


where the one-electron and two-electron integrals are 


hij — 




hijki — 


vUn)<P*iir2)^e{ri)(pkif2) 


ipj(f)dr, (15) 


Ti - r2 


dr I dr 2 ■ 


(16) 


The operators a| and aj in Eq. (14) obey the fermionic 
anti-commutation relations, 

{al,aj} = S,j, {al,a]} = {ai,aj} = 0. (17) 

In general, the Hamiltonian in Eq. (14) contains 0{N'^) 
terms, except in certain limits of very large molecules 
where use of a local baas and truncation of terms lead to 
scaling on the order of 0{N'^) [8]. The spatial encoding of 
Eq. (14) requires Q{N) qubits, one for each spin-orbital. 

While fermions are antisymmetric, indistinguishable 
particles, qubits are distinguishable and have no special 
symmetries. Accordingly, in order to construct the oper¬ 
ator SELECT (iJ), which applies terms in the second quan¬ 
tized Hamiltonian to qubits as in Eq. (2), we will need 
a mechanism for mapping the fermionic raising and low¬ 
ering operators in Eq. (14) to operators which act on 
qubits. Operators which raise or lower the state of a 
qubit are trivial to represent using Pauli matrices, 

cr+ = |1) (0| = ^ (crj - i crj) , (18) 

(J- = |0) (1| = ^ (crj -f i crj) . (19) 

Throughout this paper, trj, crj and (t| denote Pauli ma¬ 
trices acting on the j**' tensor factor. However, these 
qubit raising and lowering operators do not satisfy the 
fermionic anti-commutation relations in Eq. (17). To en¬ 
force this requirement we can apply either the Jordan- 
Wigner transformation [34, 35] or the Bravyi-Kitaev 
transformation [37-39]. 

The action of aj or Uj must also introduce a phase to 
the wavefunction which depends on the parity (i.e. sum 
modulo 2) of the occupancies of all orbitals with index 
less than j [38]. If fj G {0,1} denotes the occupancy of 
orbital j then 


«] l/Af • • • fj+1 0 fj-1 ■ ■ 

■fi) 




■ ■ fj+11 fj-1 ■ 

••/l) 

(20) 

Oj l/w • • • fj+11 fj-1 ■ ■ 

■fi) 



= (-l)^-'Mi)v 

■ ■ fj+i0fj_i ■ 

••/l) 

(21) 

a] I/a • • • fj+11 fj-1 ■ ■ 

•/i) = 0 


(22) 

Oj l/Af • • • fj+1 0 fj-1 ■ ■ 

•/i) = 0. 


(23) 


In general, two pieces of information are needed in or¬ 
der to make sure the qubit encoding of the fermionic 
state picks up the correct phase: the occupancy of the 
state and the parity of the occupancy numbers up to j. 
The Jordan-Wigner transformation maps the occupancy 
of spin-orbital j directly into the state of qubit j. Thus, 
in the Jordan-Wigner transformation, occupancy infor¬ 
mation is stored locally. However, in order to measure 
the parity of the state in this representation, one needs 
to measure the occupancies of all orbitals less than j. 
Because of this, the Jordan-Wigner transformed opera¬ 
tors are V-local, which means that some of the Jordan- 
Wigner transformed operators are tensor products of up 
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to N Pauli operators. The Jordan-Wigner transformed 
operators are 


the transformation operators in Eq. (27) and Eq. (28) so 
that 




J-l 

s — 1 

1 

" 2 


■■ ® (tI 

(24) 

i-1 

0^: = 

1 

" 2 


■ • 0 cr/. 

(25) 




I'ijkt) \q 1 q 2 q 3 q 4 ) IV') \q 1 q 2 q 3 q 4 ) ai^q^ IV') 

\ijk£) \q1q2q3q4) ak^q^ae^q^ |V') 

1-^ \ijke) \q1q2q3q4) a],q^ak,q:,ai^q^ |V’) 

1-^ \iikti \q1q2q3q4) a\,q^^],q^o-k,q 3 a(.^q^ |V') ■ ( 30 ) 


It would be convenient if we could construct 
SELECT(i7) by applying the Jordan-Wigner transform 
and acting on the quantum state, one spin-orbital in¬ 
dex at a time. For instance, SELECT(iJ) might control 
the application of a fermionic operator as follows 

\ijkt) IV') \ijk£) ai |V') 

\iik£) akQi IV') 

1-4 \ijkt) a^akai \ ip) 

1-4 \ijk£) a|a]afcaf |V’) . (26) 

However, the operators a] and aj are not unitary because 
the operators a~^ and a~ are not unitary. To correct 
this problem, we add four qubits to the selection register 
where each of the four qubits indicates whether the cr® 
or the a part of the cr^ and a~ operators should 
be applied for each of the four fermionic operators in a 
string such as ajojakai. For ease of exposition, we dehne 


new fermionic operators which are unitary, 

and aj^q 

where q € {0,1}, 



j-i 

1-1 




(27) 

S^l 



J-l 



5.0. 

o 

III 

( 8 ) 

Oj-i =icr|(^crf. 

(28) 

S = 1 




A circuit which implements these operators controlled 
on the selection register is straightforward to construct. 
Furthermore, the transformation of the terms can be ac¬ 
complished in 0(1) time. Because the Jordan-Wigner 
transformation is Wlocal, the number of gates required 
to actually apply the unitaries in SELECT(iJ) is 0{N). 
However, the terms in Eq. (27) and Eq. (28) are trivial 
to apply in parallel so that each query takes 0(1) time. 

Whereas the Jordan-Wigner transformation stores oc¬ 
cupancy information locally and parity information N- 
locally, the Bravyi-Kitaev transformation stores both 
parity and occupancy information in a number of qubits 
that scales as O(logA^) [37-39]. For this reason, the op¬ 
erators obtained using the Bravyi-Kitaev basis act on 
at most O(logiV) qubits. It might be possible to apply 
the Bravyi-Kitaev transformation with C>(logA) gates, 
which would allow for an implementation of SELECT(iJ) 
with C>(logA^) instead of 0{N) gates. However, the 
Bravyi-Kitaev transformation is much more complicated 
and this would not change the asymptotic scaling of our 
complete algorithm. The reason for this is because the 
total cost will depend on the sum of the gate count 
of SELECT(iJ) and the gate count of prepare(IT) or 
PREPARE(ri;), and the latter procedures always require 
at least 0{N) gates. 


IV. SIMULATING HAMILTONIAN 
EVOLUTION 


We use these definitions to rewrite the Hamiltonian in 
Eq. (14) so that it is explicitly a weighted sum of unitary 
Pauli products of the form we require in Eq. (1), 


^ = EE 

<?192 tj 


hij + 


+ E E 

<?192<?394 ijk£ 


02 • 


(29) 


Using the method of [33], Hamiltonian evolution can 
be simulated with an exponential improvement in preci¬ 
sion over Trotter-based methods by approximating the 
truncated Taylor series of the time evolution operator 
U = We begin by partitioning the total simula¬ 

tion time t into r segments of t/r. For each of these r 
segments we perform a Taylor expansion of the propaga¬ 
tor and truncate the series at order K, i.e. 


Inspection reveals that applying the transformations in 
Eq. (27) and Eq. (28) to Eq. (29) gives the same ex¬ 
pression as applying the transformations in Eq. (24) and 
Eq. (25) to Eq. (14). By removing factors of 1/2 from 
both transformation operators and instead placing them 
in Eq. (29), we obtain transformation operators that are 
always unitary tensor products of Pauli operators. 

Accordingly, we can implement SELECT(iJ) in the 
spirit of Eq. (26) by using four additional qubits and 


Ur = 


^ k\ 


K T / .,1 './c 

= E E k\ (31) 

fe=0 7i,'" ,7fc=l 


where in the second line we have expanded H as in 
Eq. (1). Notice that if we truncate the series at order 
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TABLE I. Database algorithm parameters and bounds 


Parameter 

Explanation 

Bound 

A 

normalization factor, Eq. (4) 

O (A^) 

r 

number of time segments, Eq. (42) 

At/ ln(2) 

K 

truncation point for Taylor series, Eq. (33) 

(n f iog(r/e) t 

(log log(r/e) J 

r 

number of terms in unitary decomposition, Eq. (1) 

O (N^) 

J 

number of ancilla qubits in selection register, Eq. (35) 

©(AlogT) 


TABLE 11. Database algorithm operators and gate counts 


Operator 

Purpose 

Gate Count 

SELECT (H) 

applies specified terms from decomposition, Eq. (2) 

0{N) 

SELECT (U) 

applies specified strings of terms, Eq. (36) 

0{NK) 

PREPARE (IT) 

prepares a superposition of states weighted by coefficients, Eq. (3) 

0{r) 

PREPARE (/3) 

prepares a superposition of states weighted by coefficients, Eq. (37) 

o{Kr) 

W 

probabilistically performs simulation under H for time t/r, Eq. (41) 

o{Kr) 

P 

projects system onto lO)®*^ state of selection register, Eq. (44) 

e(Aiogr) 

G 

amplihcation operator to implement sum of unitaries, Eq. (45) 

o{Kr) 

(PGY 

entire algorithm 

O (rKF) 


K, we incur error 


O 


{\\H\\t/rf+^ \ 
{K+l)\ )■ 


(32) 


If we wish for the total simulation to have error less than 
e, each segment must have error less than e/r. Accord¬ 
ingly, if we set r > \\H\\t then our total simulation will 
have error at most e if 




(33) 


We now discuss how one can actually implement the 
truncated evolution operator in Eq. (31). First note that 
the sum in Eq. (31) takes the form 

U = '^PjVj, j = (fc,7i,--- ,7fc), 

3 


By making 0{K) queries to SELECT(iL) from Sec¬ 
tion IV, we can implement an operator to apply the Vj 
which is referred to in [33] as SELECt(V), 

SELECT (V) \j) \lP) = \j) Vj IV') , (36) 

where the Vj are defined as in Eq. (34). This is equiv¬ 
alent to k applications of select (iL), using each of the 
| 7 i,) registers, together with k multiplications by —i. In 
order to obtain k applications of SELECT(iL), we may per¬ 
form a controlled form of SELECt(7J) K times, with each 
successive qubit in the unary representation of k as the 
control. Given that the gate count for SELECT(iL) scales 
as 0{N), we can implement SELECt(V) with 0{NK) 
gates. Applying the Pauli strings in parallel leads to 
circuit depths of 0{1) and 0{K), respectively. Table I 
lists relevant parameters along with their bounds in our 
database algorithm. Table II lists relevant operators and 
their gate counts in our database algorithm. 

We will also need an operator that we refer to as 
prepare(/3), which initializes a state. 


where the Vj are unitary and U is close to unitary. Our 
simulation uses an ancillary “selection” register |j) = 
|fc) where 0 < k < K and 1 < 7 i, < T for 

all V. We will encode k in unary, which requires Q{K) 
qubits, so that \k) = Additionally, we encode 

each | 7 „) in binary using ©(logT) qubits. While we need 
K of the | 7 t;) registers, we note that only k will actually 
be in use for a given value of \k). The total number of an- 
cilla qubits required for the selection register \j), denoted 
as J, scales as 


PREPARE (/3) jO)®'^ 


3 


(37) 


where s is a normalization factor. To implement 
prepare(/3) we first prepare the state 



{kt/rf 

k\ 


- 1/2 K 

E 





(38) 


J &Q{K\ogV)=0 


/ log [N) log (r/e) \ 
V log log (r/e) j 


(35) 


Using the convention that Ry{9) = exp[—i0cr^/2], we 
apply Ry(9i) to the first qubit of the unary encoding for 
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FIG. 1. The circuit for prepare(/3) as described in Eq. (37). 
An expression for Ok is given in Eq. (39). prepare(IF) is 
implemented using a precompnted classical database. 


k followed by Ry{9k) to the fcth qubit controlled on qubit 
k — \ for all k G [2, AT] sequentially, where 


/ 

(39) 

To each of the K remaining components of the selection 
register |7i) • • • Itr-), we apply prepare(IT) once. In 
principle, we only need to perform prepare(IT) k times, 
because the registers past k are not used. However, it 
is simpler to perform PREPARe(IT) K times, because it 
does not require control on k. 

Using results from [44], prepare(IU) can be imple¬ 
mented with (!I(r) gates by using a classically precom¬ 
puted database of the T molecular integrals. The gate 
count for prepare(/3) thus scales as 0(Kr). How¬ 
ever, this construction is naturally parallelized to depth 
0(K) + 0(r). A circuit implementing prepare(/3) is 
shown in Figure 1. 

The general strategy for implementing the truncated 
evolution operator in Eq. (34) becomes clear if we 
consider what happens to state ['(/') when we apply 
prepare(/3) followed by the operator select(U), 

SELECT (V) prepare (/3) |0)®'^|V') = V ~ ' 

(40) 

The similarity of this state to the state U |^) motivates 
the operator 

W = (prepare (/3) 0 1)^SELECT (V) (PREPARE (/3) 0 1) 

W|0)®‘'|V') = i|0)U|V') + yi-^|$), (41) 

where |<1?) is a state with the ancilla qubits orthogonal to 
jO)®'^. Note that in [33], the authors use the convention 


9k = 2 arcsin 


\ 


1 - 


{At/i 


^k-l 


(fc-l)! 


{At/rY 


that all Wiy are positive and phases are incorporated into 
the operators H^. Since we depart from that convention 
for reasons described in Section H, the second applica¬ 
tion of prepare(/3) in Eq. (41) is the transpose of the 
first application, in contrast to [33] where the conjugate 
transpose is used instead. 

At this point, we choose the number of segments to be 

r = At/\n{2). (42) 

Since A > j|i7j|, our choice of K in Eq. (33) remains 
valid. The additional factor of l/ln(2) is included to 
satisfy a requirement of oblivious amplitude amplification 
as described in [32] so that 


^ = (43) 

f k=0 

We now define a projection operator P onto the target 
state, which has support on the empty ancilla register, 

P=(jO}(Olf^0l, 

PW\Qf-^\ilj) = -\0f^U\^lj). (44) 

S 

We also define the amplification operator, 

G^-WRW'^RW, (45) 

where R = 1 — 2P is the reflection operator. With these 
definitions, we follow the procedure in [33] which uses 
the oblivious amplitude amplification procedure of [32] 
to deterministically execute the intended unitary. We 
use G in conjunction with P to amplify the target state, 

UU^li^ IV’) . (46) 

Recalling the definition of Ur in Eq. (31), our choice of 
K in Eq. (33) and r = At /In 2 imply that 

IIPG 10) 1V’) - |0)t7,lV’)|| GO (e/r), (47) 


PG10)1V’) = |0) (-U-^ 


so that the total error from applying oblivious amplitude 
amplification to all the segments will again be order e. 

The gate count of the entire algorithm is thus r 
times the cost of implementing SELECt(U) plus the cost 
of implementing PREPARe(/1). Though we implement 
SELECt(U) with 0{NK) gates, our brute-force construc¬ 
tion of prepare(IU) led to a gate count for prepare(/3) 
which scales as 0{KT). Thus, the total gate count of our 
database algorithm scales as 


O {rKT) = O 


f N^Atlog {Nt/e) \ 
V loglog(fVf/e) J 


O (N^t). (48) 


While this bound suggests an exponentially more precise 
algorithm than those based on Trotterization, in the re¬ 
mainder of our paper we discuss an even more efficient 
algorithm with improved dependence on N. 




























V. EVOLUTION UNDER INTEGRAL 
HAMILTONIANS 


In Section IV we analyzed the database algorithm for 
quantum simulating chemistry Hamiltonians in a manner 
that is exponentially more precise than Trotterization. 
The most costly part of that procedure is the implemen¬ 
tation of prepare(IV) as in Eq. (3), which prepares a 
superposition state with amplitudes that are given by 
integrals over spin-orbitals as in Eq. (15) and Eq. (16). 
Instead of classically precomputing these integrals and 
implementing prepare(IE) with a database, the strat¬ 
egy we introduce is to numerically sample the integrals 
on-the-fly using the quantum computer. Because of this, 
we call this the “on-the-fly” algorithm. To accomplish 
this, we discretize the integrals as sums and design a cir¬ 
cuit which returns the integrand of these integrals at par¬ 
ticular domain points. The motivation for approximating 
integrals as sums comes from a direct analogy between 
the discretization of time in the Taylor series approach 
for simulating time-dependent Hamiltonians [33] and the 
discretization of space in Riemann integration. 

In [33], the time-ordered exponential is approximated 
by a Taylor series up to order K, and the integrals are 
then discretized as follows on each segment. 


T exp 


nt/r 

—i H (t) dt 
Jo 


K 


k\ 


( — j-t/r 

I TH{t,)...H{h)dt 

k^O 
K 

E 

A ;=0 




k 


5] (49) 


h,---,3k=0 


Now let us suppose that our Hamiltonian does not change 
in time, but instead that the Hamiltonian itself is given 
as a definite integral over the spatial region Z so that 


H = J 'H{z)dz. (50) 

The second quantized Hamiltonian given by Eq. (29) is 
similar to this, except it includes both terms which 
are integrals over one spatial coordinate, and hijki, which 
are integrals over two spatial coordinates. While those 
integrands are also defined over all space, the integrands 
decay exponentially so we can approximate them as def¬ 
inite integrals over the finite region Z, having volume V. 
Then we can approximate the integral by 

r w 

Hzz H (z) dz zz — H {Zf,), (51) 

where Zp is a point in the domain Z at the grid point. 
As in Section IV, we begin by dividing t into r seg¬ 


ments. We turn our attention to a single segment. 



= E ^ ^ dz (52) 

k=o ■ 

where in the second line we have performed a Taylor 
expansion of the propagator and truncated at order K. 
In Eq. (52), the bolded symbol z indicates a vector of 
vectors. Like before, if r > \\H\\t then the relationship 
between K and e is given by Eq. (33). To approximate 
the integral, we divide it into p, regions of volume V/p. 
We now have 


K 

-E 


^ i-itvf '' ^ 


p^k\ 

i-itV) 


E ^ 


p=i 
fc M 


r^p^k\ 

k=0 pi,-",pfc=l 


5] H{zp,) ■■■%%,). (53) 


Eor the second quantized Hamiltonian, the in Eq. (1) 
are integrals over scalar functions w^{z) as in Eq. (5). 
Using this property it is clear that 


r 

■^ (^) = E 

7^1 

and the segment Ur can be expressed as 




i: i: 


r ^ LJj ^ A) ^ 

^ *7ir--pi,--- 

X Wj, {zp, {Zp^) ■ 




(55) 


The second quantized Hamiltonian given in Eq. (29) is 
a sum of terms which are integrals over one spatial co¬ 
ordinate and terms which are integrals over two spatial 
coordinates. This case is easily accounted for by taking 
z to be a vector including both spatial coordinates, and 
V to be the product of the volumes for the two coordi¬ 
nates. One can take the terms with the integral over the 
single spatial coordinate to also be integrated over the 
second spatial coordinate, and divided by the volume of 
integration of that second coordinate to give the correct 
normalization. We may now proceed with the truncated 
Taylor series simulation as in Section IV. Whereas our 
database algorithm required prepare(IU) to create a su¬ 
perposition of states weighted by the Wy, as in Eq. (3), 
our on-the-fly algorithm will need to create a superposi¬ 
tion of states weighted by the scalar integrands w-y{zp). 

We now show a method for dynamically decomposing 
each Wj(z) into a sum of terms which differ only by a 
sign as in Eq. (8) so that each Wj(z) is a sum of M G 
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0(inax^,y |■u;-y(i*)|/C) terms w-y^miz) G {—1,+1} where ( 
is the precision of the decomposition. First, we round 
each entry of Wy{z) to the nearest multiple of 2^ and 
decompose it in such a fashion that 


M 

'^1 (^) “ c X! 

m—1 


<C- 


(56) 


To construct prepare{w), as in Eq. (11), we perform logic 
on the output of SAMPLE(r(;) (introduced in Eq. (7), and 
detailed in Section VI) which determines whether the 
for a given value of m should be 1 or —1. Since 
the superposition in Eq. (11) must be weighted by the 
square root of this coefficient, we need to prepare states 
that either do or do not have a phase factor i so 


K) ^7 (^p)) ^ 


1 ^) ^7 (^p)) 
i \i) \w^ {Zp)) 


Wy (zp) > (2m 
w-y (zp) < (2m 


M)C 

M)C 

(57) 


where \€} = jq) \m) \p) and \w^ (zp)) was obtained from 
SAMPLe(u>). The phase factor can be obtained us¬ 
ing phase-kickback in the usual way. Then we apply 
SAMPLe(w) a second time to erase the register \wi{zp)). 
A single query to this circuit allows for the construction 
of PREPARE(?ii) with the same complexity as SAMPLe(w). 

Before explaining the integrand circuit we briefly com¬ 
ment on the additional resources required for the Taylor 
series simulation under a discretized, position-dependent, 
integrand Hamiltonian. As in the constant Hamilto¬ 
nian case, we need one register with Q{K) qubits to 
encode \k) and K registers of ©(logT) qubits to en¬ 
code |7 i)---|7 fc). However, as in the explanation be¬ 
low Eq. (11) we also need extra ancilla qubits to store 
the value of m, the grid point registers, as well as the 
value registers which are used by the integrand oracle 
SAMPLe(ui) in Eq. (7). This represents an additional an¬ 
cilla overhead of Q{K\og{Mp)). 

The sources of simulation error are also similar to the 
constant Hamiltonian case. As we show in Section VI, 
we can approximate the integrals with discrete sums to 
precision e at a cost that is logarithmic in 1/e. The error 
due to the discrete sum is controlled by the choice of 
/i, which we need to select so that the resulting error 
per segment is less than e/r. The most costly integrals 
(due to the size of their domain) will be the two-electron 
integrals in Eq. (16) which have integrands of the form 


hijk£ (^; y) 


jx) jy) (x) <pk {y) 

\x-y\ 


(58) 


where x and y represent the three spatial degrees of free¬ 
dom of two separate electrons. In Section VI, we bound 
the cost to the quantum algorithm of estimating the cor¬ 
responding integrals. 


VI. THE INTEGRAND ORACLE 

In Section V, we showed how one can implement the 
truncated Taylor series simulation technique by replacing 
a superposition state having amplitudes given by inte¬ 
grals with a superposition state having amplitudes given 
by their integrands, as well as a way of decomposing those 
integrands. We begin this section by constructing a cir¬ 
cuit which allows us to sample from the integrands as in 
Eq. (7). First, we will need a circuit which computes val¬ 
ues of the N spin-orbital basis functions (pi[zp) to piq{zp) 
at Zp, a real-valued position vector at grid point p. The 
action of each these oracles is 

Qvo \p) = Ip) {zp )), (59) 

where ^j{zp) represents the binary expansion of ‘Pj{zp) 
using logM qubits. We will need N different circuits of 
this form, one for each basis function pi{z) to (pn{z). 
Usually, the molecular spin-orbital basis functions are 
represented as sums of Gaussians multiplied by poly¬ 
nomials [43]. In that case, the circuit would be a 
reversible classical circuit that evaluates and sums the 
Gaussians associated with For example, in the 

STO-nG basis set, each orbital is a linear combination 
of n Gaussian basis functions [43]. In general, Gaussian 
functions may be evaluated classically with complexity 
that is at most polylogarithmic in 1/e [45]. The use of 
Gaussians is a historical precedent used because those 
functions are simple to integrate on a classical computer. 
However, the use of a Gaussian basis is not necessarily 
an optimal strategy for a quantum implementation. We 
leave open the problem of determining an optimal rep¬ 
resentation of molecular orbital basis functions for eval¬ 
uation on a quantum computer and develop a strategy 
based on the model chemistries used in classical treat¬ 
ments of electronic structure. 

Next, we combine N different Qy,. circuits, one for each 
y>j{z), to construct a circuit Q which allows us to apply 
any of the N basis functions. This circuit will have depth 
0(Vpolylog(Vt/e)) and may be constructed as the block 
diagonal operator 


N 

Q = W\j) {j\®Qvo- (60) 

i=i 

Thus, Q is a sequence of Qy,y circuits with the spin- 
orbital selection completely controlled on a register en¬ 
coding the basis function index j. There will a factor of 
\og{Nt/e) because the controlled operations need to ac¬ 
cess ©(logiV) qubits for j, as well as 0{\og{Nt/e)) qubits 
storing the position z. In addition, the circuit needs to 
perform analytic operations (e.g. calculating exponen¬ 
tials for STO-nG), which will contribute an additional 
factor polynomial in \og{Nt/e). 

We now discuss how one can use Q to compute the two- 
electron integrands in Eq. (58). To avoid singularities 
that would occur when two electrons occupy the same 
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point in space, we change variables in Eq. (58) so that 
^ = X — y. With this substitution, the integral becomes 

f - 0<fi{x)ipk{x - 0 ,3-,3^ 

/ - - -^- d £d X. (61) 

J 1^1 

Expressing ^ in spherical polar coordinates, with ^ = |^|, 
we have 

(fi*(x)(p*(x — ^)(pi{x)(pk{x — 0 ^ sinO d^ d9 d(j) d^x. 

(62) 

We define the maximum value of any spin-orbital func¬ 
tion as i^max and the maximum value of its derivative 
in any direction as i^max- In addition, we truncate the 
integral at a finite distance Xmax- Now assume that we 
discretize x in intervals of size 6 x along each degree of 
freedom. We can take the maximum value of ^ to be 
a^max, and choose 6 ^ = Sx, 59 = 5(j) = 5a;/a;„iax- 

The primary contribution to the complexity is in 
terms of the number of segments. The maximum 
value in the integrand of Eq. (62) is upper bounded by 
Xma.yi^fiaax- When discretizing the integral, each term in 
the sum is no larger than 0(a;niax‘/amax^®^(‘^^/^max)^) = 
^(¥’max'^a;®/a;„iax) and there are 0 {{x^ax/SxY') terms. 
Multiplying these together gives us the contribution of 
the integral to the scaling of our on-the-fiy algorithm, 

Cl (^maxa^Lx) , (63) 

which corresponds to the factor of Vmax^^, |w 3 ,(.?)| in 
Eq. (12). But how do i^max and a^max scale with N1 The 
maximum values (praax are predetermined by the model 
chemistry, and hence are independent of N. Determining 
the appropriate value of x^ax is a little more complicated. 

Because the Hamiltonian is a sum of 0{N‘^) of the inte¬ 
grals, each integral should be approximated within error 
0{e/[N'^t)) to ensure that the final error is bounded by 
e. Since the functions ^j{z) can be chosen to decay ex¬ 
ponentially, Xmax can be chosen logarithmically in the 
allowable error e. The quantum chemistry problem is al¬ 
ways defined within a particular basis, specified as part of 
a model chemistry [48]. The model chemistry prescribes 
how many spin-orbitals, how many basis functions, and 
what type of basis functions should be associated with 
each atom in a molecule. This includes a specification for 
parameters of the basis functions which impose a particu¬ 
lar maximum value (/^max) as well as a cutoff distance be¬ 
yond which each is negligibly small. However, the 

space between basis functions on different atoms must 
grow as the cube root of N, because the molecular vol¬ 
ume will grow as 0{N). This would imply that the value 
of Xmax needed scales as 

a^max G O log (Nt/e)^ . (64) 

Nevertheless, each individual orbital <fj{z) is non- 
negligible on a region that grows only as O(logiV) for 


a given model chemistry. It is therefore advantageous 
to modify the grid used for the integral so it only in¬ 
cludes points where one of the associated orbitals is non- 
negligible. This can be performed at unit cost if the 
center of each spin-orbital function is provided in an ad¬ 
ditional register when querying the circuit Q. As above, 
the region where the orbital can be regarded as non- 
negligible can be chosen logarithmically in e, to ensure 
that the overall error in the simulation is within error e. 

To be more specific, the coordinates for x can be cho¬ 
sen to be centered around the center of orbital ipi, with 
the components of x only going up to a maximum value 
scaling as 

Xmax G O (log (Nt/e)). (65) 

Eor we only wish to take values such that ipj{x — ^) are 
non-negligible. Here it should be noted that the spheri¬ 
cal polar coordinates are only advantageous if we are in 
a region where ^ is near zero, where the Cartesian coor¬ 
dinates would have a divergence. In regions where ^ is 
large, the extra factor of ^ for the integral in spherical 
polar coordinates increases the complexity. 

Therefore, if the minimum value of j^j such that 
(pj{x — ^) is non-negligible is O (log (Nt/e)), then the 
maximum value of j^j such that ipj{x—^) is non-negligible 
will also be O (log (Nt/e)). Therefore we can use spher¬ 
ical polar coordinates, and obtain scaling as in Eq. (63) 
with Xmax as in Eq. (65). On the other hand, if the min¬ 
imum value of 1^1 such that (pj{x — ^) is non-negligible 
is O (log (Nt/e)), then we can use Cartesian coordinates, 
and the division by j^j can only lower the complexity. 
We obtain a contribution to the complexity scaling as 
^ (‘t^max^max) with Xmax as in Eq. (65). Here the power 
of Xmax is 3 rather than 5, because we divide instead of 
multiplying by j^j as we did spherical polar coordinates. 

Next we consider the grid size needed to appropriately 
bound the error in the discretized integration. The anal¬ 
ysis in the case where Cartesian coordinates are used is 
relatively straightforward. Considering a single block in 
6 dimensions with sides of length Sx, the value of the 
integrand can only vary by the maximum derivative of 
the integrand times 6 x (up to a constant factor). The er¬ 
ror for the approximation of the integral on this cube is 
therefore that maximum derivative times <5x^. Then the 
number of these blocks in the integral is C>((xmax/<Ix)®), 
giving an overall error scaling as x%,^y.5x times the max¬ 
imum derivative of the integrand. 

The maximum derivative of the integrand can be 
bounded in the following way. For the derivative with 
respect to any component of x, we obtain the derivative 
of the integrand scaling as 

O (^ (66) 

\ (^niax / 

where we have used the fact that we are only using Carte- 
sian coordinates for j^j = D(xmax)- For the derivative of 
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the integrand with respect to any component of ^ in the 
numerator of the integrand, the same scaling is obtained. 
For derivatives with respect to components of ^ in the de¬ 
nominator, the scaling is 

O (. (67) 

Va^max/ 

Overall, we therefore bound the error when discretizing 
in Cartesian coordinates as 

^ ((V^max + V^max/a^max) <7’maxa^max'^a;) . (68) 

The analysis for spherical polar coordinates is a little 
more subtle, but it is largely equivalent if we scale the an¬ 
gular variables. It is convenient to define scaled angular 
variables 


h — a^max^ 


(69) 


Then the discretization lengths for all variables are 5x. 
The volume of each block in the discretization is again 
and there are 0((a:niax/^a:)®) blocks. The total er¬ 
ror is again therefore the maximum derivative of the in¬ 
tegrand multiplied by x%^^y.6x. 

The derivative of the integrand with respect to any 
component of x is again given by Eq. (66). Multipli¬ 
cation by ^ gives a factor ©(Xmax), but the change of 
variables to 6 ' and (j)' gives division by a factor of s^ax- 
The derivative of the integrand with respect to O' or 0' 
in any of the spin orbitals again gives a factor scaling as 
in Eq. (66) . The derivative of the integrand with respect 
to ^ or O' in ^sin(0'/a:max) scales as in Eq. (67). 

As a result, regardless of whether Cartesian coordi¬ 
nates are used or spherical polar coordinates, the error 
due to discretization is bounded as in (68). Thus, to 
achieve error in the integral no larger than 0 {e/{NH))^ 
we require that 


5x G Q 


(T^max “f V^max/a^niax) 'T’max^^max 


(70) 


The total number of terms in the sum then scales as 



= 0 


' /NH 


(‘T’max + T^max/ ^max)^: 


X^ 

maxmax 


(71) 


This is quite large, but because we only need to use a 
number of qubits that is the logarithm of this, it only 
contributes a logarithmic factor to the complexity. Be¬ 
cause the logarithm scales as C>(log(A^t/e)), it contributes 
this factor to the complexity of SAMPLE(r(;). 

Given Q, computing the integrand in Eq. (62) is 
straightforward. We need to call Q four times on reg¬ 
isters that contain x and f. We then multiply those out¬ 
puts together with the value of ^ sin0. Next we consider 
how to construct a circuit for the one-electron integrals 


in Eq. (15). First, one constructs N additional circuits 
similar to the ones in Eq. (59) that return W‘^(pj{z) as 
opposed to ^j{z). These oracles are incorporated into a 
one-electron version of Q which is called along with a rou¬ 
tine to compute the nuclear Coulomb interactions. The 
one-electron integrals have singularities at the positions 
of the nuclei. Similar to the two-electron integrals, these 
singularities can be avoided by using spherical polar co¬ 
ordinates. Each term in the sum over the nuclei should 
use spherical polar coordinates centered at that nucleus. 
Selection between the one-electron and two-electron rou¬ 
tines is specified by jq). Putting this together, we can im¬ 
plement SAMPLE(iy) as in Eq. (7) with C>(A^log A^) gates, 
and, as discussed in Section V, prepare('u;) has the same 
complexity. 

While the 0{N) gate count of prepare('u;) is much less 
than the 0 {N'^) gate count of prepare(IT), our on-the- 
fiy algorithm requires more segments than the database 
algorithm. Whereas our database algorithm requires 
r = At/ln(2) segments where A is the normalization in 
Eq. (3), our on-the-fly algorithm requires r = At/ln(2) 
segments where A € ©(Ft/s^ax^max) is the normaliza¬ 
tion in Eq. (11), which is accounted for in Eq. (12) and 
Eq. (63). Thus, by performing the algorithm in Sec¬ 
tion IV using prepare(i(;) instead of prepare(IT) and 
taking r = A</ln(2), we see that our on-the-fly algorithm 
scales as 


O (rNK) = O (NKXt). (72) 

Using the scaling in Eq. (65), we can bound A as 

A e O (r(^Lx^Lx) e o (iV^log iNt/e)f) . (73) 

so that the overall gate count of the on-the-fiy algorithm 
scales as 


O (N^Kt) = O (N^t) . (74) 

Recall that the O notation indicates that logarithmic fac¬ 
tors have been omitted. The full scaling includes a power 
of the logarithm of 1/e. 


VII. DISCUSSION 

We have introduced two novel algorithms for the sim¬ 
ulation of molecular systems based primarily on the re¬ 
sults of [33]. Our database algorithm involves using a 
database to store the molecular integrals; its gate count 
scales as 0{N^t). Our on-the-fiy algorithm involves com¬ 
puting those integrals on-the-fly; its gate count scales 
as 0{N^t). Both represent an exponential improvement 
in precision over Trotter-based methods which scale as 
0 {N^y/t^/e) when using reasonably low-order decompo¬ 
sitions, and over all other approaches to date. 

Specifically, our database algorithm scales like 
0{N^At) where we have used the bound A S 0{N'^). 
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However, we believe this bound is very loose. As dis¬ 
cussed in [8, 43], the use of local basis sets leads to a 
number of two-electron integrals that scales as 0 {N‘^) in 
certain limits of large molecules. Accordingly, the true 
scaling of the database algorithm is likely to be closer 
to 0{N^t). It also seems possible that our integration 
scheme is suboptimal; it is possible that it can be im¬ 
proved by taking account of smaller values of 

Our asymptotic analysis suggests that these algorithms 
will allow for the quantum simulation of molecular sys¬ 
tems larger than would be possible using Trotter-based 
methods. However, numerical simulations will be crucial 
in order to further optimize these algorithms and bet¬ 
ter understand their scaling properties. Just as recent 
work showed significantly more efficient implementations 
of the original Trotterized quantum chemistry algorithm 
[5-9], we believe the implementations discussed here are 


far from optimal. Furthermore, just as was observed for 
Trotterized quantum chemistry [7, 9] , we believe our sim¬ 
ulations might scale much better when only trying to 
simulate ground states of real molecules. In light of this, 
numerical simulations may indicate that the scaling for 
real molecules is much less than our bounds predict. 
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